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Abstract 

We present the first results of time dependent two-fluid cosmic-ray (CR) modified MHD 
shock simulations. The calculations were carried out with a new numerical code for 1-D 
ideal MHD. By coupling this code with the CR energy transport equation we can simulate 
the time-dependent evolution of MHD shocks including the acceleration of the CR and 
their feedback on the shock structures. We report tests of the combined numerical method 
including comparisons with analytical steady state results published earlier by Webb, as 
well as internal consistency checks for more general MHD CR shock structures after they 
appear to have converged to dynamical steady states. We also present results from an 
initial time dependent simulation which extends the parameter space domain of previous 
analytical models. These new results support Webb's suggestion that equilibrium oblique 
shocks are less effective than parallel shocks in the acceleration of CR. However, for realistic 
models of anisotropic C’R diffusion, oblique shocks may achieve dynamical equilibrium on 
shorter timescales than parallel shocks. 
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1. Introduction 


Nonlinear theories of diffusive shock acceleration have demonstrated the importance of 
cosmic-rav (CR) feedback on the evolution of shock structures (e.g., Blandford k Eichler 
1987). Using two-fluid models it has been shown (e.g. Drury& Volk 1981, Achterberg 
e.t al. 1984. Kang k Jones 1990) that CR pressures can become large enough to 
smooth shocks, eliminating the entropy generating gas sub-shock. More generally, the CR 
feedback modifies the efficiency' of energy transfer from gas to CR in the shock. Recent 
numerical simulations (e.g., Drury 1' Falle 1987, Falle V' Giddings 1989, Jones & Kang 
1990. Kang k Jones 1991, Kang. Jones k Ryu 1992) have also shown the importance 
of time-dependent effects in the determination of CR modified shock properties. For the 
most part, past studies of CR modified shocks have focused on pure hydrodynamical, 
(or parallel, sonic mode) models of the shock dynamics. Magnetic fields, however, will 
generally lie dynamically important in many environments where particle acceleration 
occurs. It has been suggested that, components of the magnetic field which axe aligned 
perpendicular to the shock normal (tangent to the shock face) can alter the efficiency 
of the acceleration process. Jokipii (1987) has pointed out that, for standard models 
of CR anisotropic diffusion (see equation [4.1] below) perpendicular components of the 
field will decrease the shock crossing time for a CR particle, increasing the rate at which 
individual particles gain energy from the shock. On the other hand, Baring, Ellison k Jones 
(1993) have shown that the efficiency of thermal particle injection into the CR population 
can be dramatically decreased by perpendicular magnetic field components, i From these 
examples it is clear that to understand tin* fundamental nature of shock acceleration in 
more realistic astrophysical settings, full MHD calculations are needed. Two-fluid models 
of C’R transport along the lines introduced by Drury V Volk ( 1981 ) are an efficient means to 
begin such explorations. Such models enable one to economically calculate the dynamical 
features of flows within the constraints imposed by the need to estimate a prion some 
closure parameters for the CR. Equilibrium MHD CR-modified shock structures have been 
calculated bv Webb (1983) using these methods for the case where the gas is cold and its 
pressure call be ignored. Webb's calculations demonstrated that, as for the CR-modified 
gas dynamic flows, both sub-shock and smoothed. CR dominated solutions to the MHD CR 
shock equations were possible. However, his solutions suggested that, m the limited range 
of conditions he could consider, shock acceleration of cosmic-rays was less effective when 
the upstream tangential components of the field were strong. Among other consequences 
this appears to expand the portions of the shock parameter space that lead to sub-shock 
solutions. That could impact on such issues as low energy injection processes and the 
momentum distribution of the CR. 


Some subsequent steady state two-fluid analyses, considering a wider parameter space 
support the above impressions ( Kennel e.t al 1985. Webb et al 1986). In this paper we report 
the first results of time dependent MHD CR two-fluid simulations. We present tests o 
a new numerical code against Webb's analytical models as well as more general internal 
checks on the code’s ability to evolve shocks to self-consistent, steady state MHD CR 
structures. Our simulations confirm Webb's calculations. Our numerical models also allow 
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us to extend Webb's calculations by lifting the cold gas. (P, = 0). restriction to examine 
the effects of a full range of initial parameters on CR-modified MHD shock structures. 

2. Methods 


We solve the equations of ideal MHD for one dimensional flow in Cartesian coordinates 
(e.g.. Jeffrey 196S). As with two-fluid gas dynamic models the conservation equations are 
modified to include momentum and energy source terms from CR feedback (e.g., rury ' 
Volk, 1981. Jones &: Kang. 1990). The MHD equations are written in conservative, vector 


form as 


with 


h> + ^I = s . 

dt dx 


l p \ 

(»'X 

P ll y 
q= | P«z 

U IJ 

B- 

V E ) 


( 2 . 1 ) 


(2.2) 


and 


/ 


pu.r 


\ 


pi'i + p* - B‘x 

pa x tt y — B.i-By 
pii.rttz — BxBz 
By ux — Bjtty 

. B’zttx - Bxu : , 

\(E + p* ) It X - B.r(Bxtl.V + ByUy + B=U Z )J 


(2.3 


while the CR source term vector is 
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The CR are themselves treated as a massless, diffusive fluid through a conservation 
equation for the CR energy. £,. derived from the diffusion-advection equation (Skilling 

1975): namely. 
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In these relations 
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Pc = he - 1 )Ec. 

and the magnetic field components are expressed in rationalized units 
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In the expressions presented above the following definitions hold: p is the mass density; 
I( B x and a,,. By, B : are the components of velocity and magnetic field parallel 
and perpendicular to shock front: P,,, ~ iy and P c , ~,c are the gas and cosmic-ray pressures 
and adiabatic indices, and k is the energy weighted spatial diffusion coefficient for the CK 
parallel to the shock normal. The quantity S f is a. term that allows direct energy transfer 
from the °as to the CR. such as through the low energy injection of thermal particles into 
the CR population (e.g.. Jones & Kang 1990). This is introduced for completeness, but 
for the present, we set S e = 0. In this discussion we set lg = 5/3, while 7c , which depends 
upon the mix of nonrelativist ic and relativistic particles in the CR population, will be 
treated as an input parameter. In general both -, c and k are properties of the solution so 
the need to specify their properties a priori, is the major drawback of the two-fluid model. 


In the interest of simplicity we will not consider here flows with tangential magnetic 
field rotations, although the numerical code is quite capable of handling such features. 
Thus, without further loss of generality we can place the magnetic field in the X-Z plane 
jj _ ( £? r . 0. ). We will also restrict ourselves to flows with no upstream tangential 

velocity. *u u = u : = 0- All of the simulations discussed here are essentially piston driven 
shock tubes. \Ye establish the flows by projecting magnetized fluid with embedded CR 


m from the right boundary, using an open boundary condition and reflecting it off a 
wall (piston) at the left boundary. The tangential magnetic field at the left boundary is 
"mirrored” in the same manner as the gas density. Previous simulations of CR modified 
shocks have shown that their time to evolve to dynamical equilibrium scales with the so- 
called diffusion time. t d . (e.g. Jones X Kang 1990). which in the present case is conveniently 

expressed as 
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where a, is the shock speed (see equation. 3.2). In our discussion below we will express 
simulation times in units of t (h appropriate to that simulation. The width of a CR-modified 
shock transition scales with the related diffusion length scale. x d . 
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To obtain accurate numerical results with the methods we employ it is important that 
computation zones be small enough to resolve the diffusive shock features on this scale. 
For our discussion we define, therefore, the resolution ratio of each simulation to be, 


Our numerical method solves equation [2.1] through the second order finite difference 
-Total Variation Diminishing’- (TVD) gas dynamics method of Harten (Harten 1981) 



extended to MHD (Ryu k Jones 1993). The pure MHD (S = 0) form of the equation 
is solve with the aid of an approximate Rieinann solver, used to estimate the fluxes, F. 
CR source corrections, S, are then added in a manner preserving second order accuracy. 
Shocks and other discontinuities are generally resolved within a few zones. The CR energy 
equation is solved using a second order combined monotone advection and Crank-Nicholson 
scheme. Further details of the method will be presented elsewhere (Frank, Jones k Ryu, in 
preparation). A pure gas dynamical version of the code was tested against both analytical 
steady state solutions and numerical time dependent models calculated with our well- 
tested PPM code (Jones k Kang 1990). with excellent agreement. The pure MHD code 
was tested against a variety of 1-D shock tube problems involving all three families of 
MHD waves vising a nonlinear MHD Rieinann solver (Ryu k Jones 1993). 

3. Results: Comparisons with Analytical Models 


In order to test the accuracy of our numerical method we first attempted to reproduce 
the analytical steady state solutions of Webb (19S3). In that paper Webb demonstrated 
in addition to discontinuous gas sub-shock solutions with a smooth CR shock precursor, 
that one may obtain completely smooth shock solutions to the MHD CR equations if 
the downstream velocity remains super-Alfvenic and the upstream CR pressure, P c , is 
high. That behavior is analogous to smooth gas dynamic shock solutions identified earner 
by Drury t Volk (19S1). However. Webb was able to consider only flows in which the 
upstream gas was cold: i .<■*.. in which P g = 0. 

In figures 1 and 2 wo present the rime-asymptotic shock structures formed in numerical 
simulations with upstream conditions corresponding to those m Webb s paper (his figure < 
and figure S). The upstream conditions for these simulations are given, as models 1 and 2, in 
table 1 . These simulations were performed with a constant CR diffusion coefficient, k = .01. 
The resolution ratios. /?,•• for models 1 and 2 are n r = 21 and 32 respectively. The Alfvenic 
Mach numbers for models 1 and 2 are M a = u p} /p/B x - 1 and 2 respectively where u p is 
the piston speed. In these models = 4/3. The simulations where carried out until the 
postshock state appeared steady; namely, t ~ 30 tj. The resultant shock transformations 
provide excellent agreement with our best estimates of the downstream states found b\ 
Webb. The largest uncertainties in the comparison are. in fact, the detei urination of the 
downstream states from Webb's figures. For model 2 the entropy, s = log(P g / p 1 * ), (not 
shown in the figure) increases though the shock, demonstrating that the flow contains an 
MHD fast mode sub-shock, as predicted by Webb. The entropy for the flow in figure 1 
shows no increase, again as predicted. Recall that these particular calculations were meant 
to be carried out under Webb's cold gas (P g « 0) restriction. For numerical reasons the 
upstream gas pressure was set in practice to be a small fraction (10 **) of the CR pressure. 
As an additional comparison we have also reproduced the pure perpendicular (B x — 0) 
smooth and MHD sub-shock models of Webb with the accuracy comparable to that shown 
in figures 1 and 2. 

Our previous simulations of CR-modified gas shocks required numerical grids that over 
resolved the CR shock precursor by roughly a factor of 10 to assure high accuracy (Jones & 
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Kang 1990). In order to explore the dependence of the accuracy on numerical resolution for 
these MHD simulations we have run a series of tests with the upstream conditions of model 
2. varying the numerical resolution. We ran simulations with n r = 4,8,16,32. In figure 3 
we plot a measure of the fractional error in the downstream CR pressure, (compared with 
the value for the highest resolution case. ??,. = 32), 


Cc = P c (n,.)-P c ( 32) 
Pc( 32). 


(3.1) 


The Figure shows that the simulated shocks converge quickly once n r > 10. As mentioned 
above, that state is in good agreement with Webb’s analytical result. The converged 
numerical postshock CR pressure that is about d'/c higher than oui estimate of Webb s 
result, although we attribute much of this error to uncertainties in reading final states from 
published figures rather than tables. These results agree well with previously mentioned 
gas dynamic behaviors found by Jones k Kang (1990). We believe the limiting factor 
to be the accuracy of the Crank-Nicholson method used for updating the diffusive CR 
energv equation. As a further test of the numerical method we have perfoimed tests of the 
self- consist enc v of more general steady state klHD CR shock solutions. This was done by 
testing the accuracy of various klHD jump conditions for apparently steady shocks in the 
shock frame. For example, the various momentum components of the full flux vector. F, 
in equation [2.3] should be the same across the shock when measured in that frame (with 
F '2 corrected to include P ( ) The shock velocity for this exercise was determined from 
the conservation of mass equation for a steady shock transformed into the piston frame, 
namely. 

= M (3.2) 

[PI 

where [ j refers to differences across tin' shock. V\e find that for the simulations with 
?/ r > 10 the jump conditions were satisfied to better than one part in 1 x 10 . 


4 Finite Gas Temperature MHD Shocks 


Since our numerical code solves the full MHD CR two-fluid equations there is no need 
to restrict investigations to those cases where P ( j — 0. In figure 4 we piesent the results 
from a simulation of an MHD CR shock formed from gas of finite upstream pressure, P g . 
The upstream conditions for this simulation, model 4, is presented in table 1. We note 
that the sonic Mach number in this case is M s = 4. The Alfvenic Mach number M a = 12. 
Standard weak scattering models of particle diffusion lead to differences in diffusion across 
and along field lines. Thus k, which controls diffusion along the shock noimal should 
depend on the angle between field and the shock normal, d = tan Thus we adopt 

k of a form, (e.g., Webb 1983, Jokipii 198*, Zank et. al. 1990), 

s — k || cos"* o T k j_ sin' o, ( ^ ) 

where the directions || and _L refer to the magnetic field direction. Note, of course, that <j> is 
generally a changing function of space and time. For this initial exploration we arbitrarily 
set k || = .01 — 10 x 
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In general P c develops more rapidly for a stiffer CR equation of state (Jones &: Kang 
1990). Thus, in order to keep the computational costs low for these first tests we used 
~ IC = 5/3 in the following model. This would influence the detailed properties of the steady 
state flow, hut will not alter the qualitative character in ways that are important to the 
present discussion. 

In figure 4 we illust rate the evolution of a time dependent, finite gas pressure simulation. 
In this model the upstream magnetic field angle is o 0 = 30°. We present results at three 
different times: t = 12,24. and 3 6f f /. Since k is not a constant in space or time in this 
simulation, we define t d here in terms of k = 0.01. Figure 4 shows a strong fast mode 
shock driven by the piston. That shock has become almost smoothed by the CR pressure. 
Comparisons of this model with an analogous parallel field (<p = 0°) simulations show a 
number important differences. First, the parallel shock model reaches a dynamical steady 
state more slowly than the oblique shock case. That is simply because, according to 
equation [4.1] the diffusion coefficient, k. in the parallel case is greater, so that energy 
gain by the C'R is slowed (e.g.. Jokipii 19S7). On the other hand, while the dynamical 
steady state may be reached more quickly in oblique shocks than in parallel shocks the 
effectiveness of the acceleration in the oblique case is reduced. The downstream value of 
P r in the oblique shock case is decreased by S% from what is obtained in the parallel shock 
model. It is reasonable to expect that in the oblique models the upstream momentum flux 
which would have gone into accelerating CR is being used to do work on the tangential 
magnetic field. We note that behind the fast mode shock a weaker slow mode shock 
compresses the CR driven transition density spike (see Jones <Sc Kang 1990 for a discussion 
of this feature)- Close examination also shows that CR particle acceleration is taking place 
as that slow mode shock develops. Modification of the density spike, which can only be 
seen in time dependent models, is an example the additional complications which arise due 
to the multiple wave families present, in MHD. 

5. Conclusions 


1) The numerical code we have developed accurately computes two-fluid models of 
CR-modified MHD shocks. If the resolution ratio defined in the text, n r > 10, then the 
time asymptotic properties of the simulations appear to converge to analytically predicted 
steady states. The time asvmptotic numerical shocks are also internally consistent in terms 
of conservation laws expected to lx* satisfied across steady shocks to at least one part in 
10 3 4 . 


2) Because of the work done on tangential fields within the shocks, time asymptotic 
particle acceleration will tend to be more efficient in parallel shocks than in oblique 
shocks. However, the oblique shocks reach dynamical steady states more quickly for a 
given diffusion coefficient parallel to the magnetic field. 

3) Transient features develop from the MHD fast mode CR shocks similar to those seen 
in pure CR hydrodynamical shocks. However, these transients are modified and made more 

complex by the development of MHD slow mode shocks. 
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Figure Captious 


Pig i. — Model 1 in table 1. MHD CR shock transition region for a piston driven shock 
with upstream conditions taken from figure < of Webb (19S3). Shown are the density, 
p. normal component of velocity, n x . tangential velocity, u :> tangential component of 
magnetic field. B : . magnetic field orientation angle, o = ton and cosnnc-ray 

pressure. P c . See table 1 for upstream flow conditions. The abscissa is given in units of 
diffusion length x (t = k/m*. The dashed lines are post-shock values taken from Webb's 
figure 7. except in the plot of u : where the value was calculated using equation [3.2] 
and the MHD steady state jump conditions. 

Fig. 2. Model 2 from table 1. MHD CR Shock transition region for a piston driven 
shock with upstream conditions taken from Webb s figure S. 

Fig. 3. Fractional error (equation [3.1] in the computed value of the post-shock CR 
pressure. P r . in model 1 as a function of resolution, n r = .r^/A.r. 

Fig. 4.- Model 3 from table 1. This model has an upstream field orientation, <f> 0 = 30°. 
Results are shown at t — 12.24 and 3 Gt ( /. 
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